library(here)
library(haven)
library(ggplot2)
library(dplyr)
library(plm)
library(lmtest)
library(magrittr)
library(raster)
library(gridExtra)
library(stargazer)
library(ggridges)
library(grid)


# Data -------------------------------------------------------------------

# Main data
df1 <- read_dta(here("Data", "adm3_panel.dta"))
df1$time <- df1$year - 1991

# Collapsed by fragmentation type
dfc <- read_dta(here("Data", "bfcollapse.dta"))


# Descriptives ------------------------------------------------------------

# Lights over time, by fragmentation type
dfc <- dfc %>% group_by(year) %>% mutate(non = max(lights))

# Figure 4
ggplot(dfc, aes(y=lights, x=year)) + 
  geom_line(data=dfc, size=1.2, 
            aes(linetype=factor(split_type, labels=c("None","Mother","Splinter")),
                color=factor(split_type))) +
  labs(x=" ", y="Lights", linetype="Fragmentation Type") +
  scale_x_continuous(limits=c(1992,2013), breaks = round(seq(min(dfc$year), max(dfc$year), by = 3),1)) +
  scale_color_manual(values=c("gray15","gray15","gray15"), guide="none") +
  scale_linetype_manual(values=c(1,6,3)) +
  guides(linetype = guide_legend(override.aes=list(color=c("gray15","gray15","gray15")))) +
  geom_vline(xintercept=1996, color="black",alpha=.5,size=1) +
  geom_vline(xintercept=2001, color="black",alpha=.5,size=1) +
  theme_bw() +
  theme(panel.grid = element_blank(),
        legend.position = "bottom") +
  annotate("text", x = 1994.5, y = .35, label = "Fragmentation (1996)", size = 2) +
  annotate("text", x = 2004.3, y = .35, label = "Equal periods pre- and post-fragmentation (2001)", size = 2)
ggsave(file = "bftsplot.png", path = here("Figures"), device = "png", width = 7, height = 3.5, units = "in", dpi = 500)


# Lights density plots (Figure A3)
ggplot(df1, aes(lights)) + 
  geom_density(alpha=.5, colour="black", fill = "gray85") + 
  labs(x="Lights", y="Density") + 
  theme_bw() + 
  theme(panel.grid.minor = element_blank())
ggsave(file="bfdens1.png", path = here("Figures") ,device = "png", width=3.5, height=3.5,units="in", dpi = 500)

ggplot(df1, aes(ln_lights)) + 
  geom_density(fill="gray85", alpha=.5, colour="black") + 
  labs(x="log(Lights + 0.001)", y="Density") + 
  theme_bw() +
  theme(panel.grid.minor = element_blank())
ggsave(file="bfdens2.png", path = here("Figures") ,device = "png", width=3.5, height=3.5,units="in", dpi = 500)


# Summary table (Table A1)
dfSum <- dplyr::select(df1, lights, ln_lights, motherpost97, splinterpost97, l_ln_pop, l_spei)

stargazer::stargazer(as.data.frame(dfSum), 
                     covariate.labels = c("Lights", "log(Lights)","Mother","Splinter","log(Population Density)","Drought Index"),
                     title="Summary Statistics",
                     label="sumStats",
                     out=here("Tables", "sumStats.tex"),
                     notes = c("Note: sample includes 351 departments."))




# Diff-in-Diff ------------------------------------------------------------

# Create empty vector to hold results
mainMods <- vector("list", 3)

# All models are estimated using the PLM package
# model="within" for within department effects
# effect="twoays" for both departement and year FEs
# index=c("id_3","year") to define panel


##### diff diff, no controls #####
mainMods[[1]] <- plm(ln_lights ~ motherpost97 + splinterpost97,
                     data = df1,
                     model = "within",
                     effect="twoways",
                     index=c("id_3", "year"))

##### diff diff, pop spei #####
mainMods[[2]] <- plm(ln_lights ~ motherpost97 + splinterpost97 + l_ln_pop + l_spei,
                     data = df1,
                     model = "within",
                     effect="twoways",
                     index=c("id_3", "year"))

##### diff diff, pop spei trend #####
mainMods[[3]] <- plm(ln_lights ~ motherpost97 + splinterpost97 + time*(l_ln_pop + l_spei),
                     data = df1,
                     model = "within",
                     effect="twoways",
                     index=c("id_3", "year"))


# Empty vectors to capture elements for clustered standard errors correction 
# Clustered HC1 standard errors
mainCladj <- vector("list",3)
cvcov <- vector("list", 3)
clMain <- vector("list", 3)

for(i in 1:3){
  
  G <- length(unique(df1$id_2))
  N <- length(df1$id_2)
  mainCladj[[i]] <- (G/(G-1)*(N-1))/mainMods[[i]]$df.residual
  cvcov[[i]] <- mainCladj[[i]]*vcovHC(mainMods[[i]], type="HC1", cluster="group", adjust=T)
  clMain[[i]] <- coeftest(mainMods[[i]], vcov = cvcov[[i]])
  
}


# Main coefficient plot 
mainPlot <- data.frame(rbind(
  cbind(clMain[[1]][1:2], clMain[[1]][1:2,2], "3No controls"),
  cbind(clMain[[2]][1:2], clMain[[2]][1:2,2], "2Controls"), 
  cbind(clMain[[3]][1:2], clMain[[3]][1:2,2], "1Controls-Trend")))
mainPlot$Var <- c("Mother", "Splinter")
colnames(mainPlot) <- c("Estimate", "SE", "Model", "Var")
mainPlot[,1:2] %<>% lapply(function(x) as.numeric(as.character(x))) # Convert several columns to numeric from factor

ggplot(mainPlot, aes(y=Estimate, x=Model)) +
  geom_linerange(aes(ymin=Estimate - 1.96*SE, ymax=Estimate + 1.96*SE), lwd=1) +
  geom_pointrange(aes(ymin=Estimate - 1.68*SE, ymax=Estimate+1.68*SE), lwd=1.5, fill="white", shape=21) +
  coord_flip() + facet_wrap(~Var,ncol=1) + labs(x="") + 
  theme_minimal() +
  geom_hline(yintercept=0, size=1.2, color="darkred", alpha=.5) + 
  scale_y_continuous(limits=c(-0.75, 0.2)) +
  scale_x_discrete(labels=c("Controls-Trend", "Controls", "No Controls")) +
  theme(strip.text = element_text(size=11),
        panel.grid.minor.x  = element_blank(),
        panel.grid.major.x  = element_blank())
ggsave(file = "mainModsplot.png", path = here("Figures"), device = "png", width = 4.5, height = 3.5, units = "in", dpi = 500, cplot1)



##############################
##### Subset to pre 2002 #####
##############################

df2 <- filter(df1, year < 2002)

# Empty vector to store results
subMods <- vector("list", 3)

##### diff diff, no controls #####
subMods[[1]] <- plm(ln_lights ~ motherpost97 + splinterpost97,
                    data = df2,
                    model = "within",
                    effect="twoways",
                    index=c("id_3", "year"))

##### diff diff, pop spei #####
subMods[[2]] <- plm(ln_lights ~  motherpost97 + splinterpost97 + l_ln_pop + l_spei,
                    data = df2, 
                    model = "within",
                    effect="twoways",
                    index=c("id_3", "year"))

##### diff diff, pop spei trend #####
subMods[[3]] <- plm(ln_lights ~  motherpost97 + splinterpost97 + time*(l_ln_pop + l_spei),
                    data = df2,
                    model = "within",
                    effect="twoways",
                    index=c("id_3", "year"))

# Cluster SEs
subCladj <- vector("list",3)
subcvcov <- vector("list", 3)
clSub <- vector("list", 3)

for(i in 1:3){
  
  G <- length(unique(df2$id_2))
  N <- length(df2$id_2)
  subCladj[[i]] <- (G/(G-1)*(N-1))/subMods[[i]]$df.residual
  subcvcov[[i]] <- subCladj[[i]]*vcovHC(subMods[[i]], type="HC1", cluster="group", adjust=T)
  clSub[[i]] <- coeftest(subMods[[i]], vcov = subcvcov[[i]])
  
}


subPlot <- data.frame(rbind(
  cbind(clSub[[1]][1:2], clSub[[1]][1:2,2], "3No controls"),
  cbind(clSub[[2]][1:2], clSub[[2]][1:2,2], "2Controls"), 
  cbind(clSub[[3]][1:2], clSub[[3]][1:2,2], "1Controls-Trend")))
subPlot$Var <- c("Mother", "Splinter")
colnames(subPlot) <- c("Estimate", "SE", "Model", "Var")
subPlot[,1:2] %<>% lapply(function(x) as.numeric(as.character(x)))


ggplot(subPlot, aes(y=Estimate, x=Model)) +
  geom_linerange(aes(ymin=Estimate - 1.96*SE, ymax=Estimate + 1.96*SE), lwd=1) +
  geom_pointrange(aes(ymin=Estimate - 1.68*SE, ymax=Estimate+1.68*SE), lwd=1.5, fill="white", shape=21) +
  coord_flip() + facet_wrap(~Var,ncol=1) + labs(x="") + 
  theme_minimal() +
  geom_hline(yintercept=0, size=1.2, color="darkred", alpha=.5) + 
  scale_y_continuous(limits=c(-0.75, 0.2)) +
  scale_x_discrete(labels=c("Controls-Trend", "Controls", "No Controls")) +
  theme(strip.text = element_text(size=11),
        panel.grid.minor.x  = element_blank(),
        panel.grid.major.x  = element_blank())
ggsave(file = "subModsplot.png", path=here("Figures"), device = "png", width=4.5, height=3.5, units = "in", dpi = 500)





# Main regression table (Table A2)

stargazer(mainMods[1], mainMods[2], mainMods[3], subMods[1],subMods[2],subMods[3], 
          se = list(clMain[1][[1]][,2], clMain[2][[1]][,2], clMain[3][[1]][,2], clSub[1][[1]][,2],clSub[2][[1]][,2],clSub[3][[1]][,2]),
          omit.stat = c("f","adj.rsq"),
          dep.var.labels = "Lights",
          label = "bigtable",
          title="Estimated effect of province fragmentation on department light intensity",
          covariate.labels = c("Mother","Splinter","Population","Drought Index","Trend:Population","Trend:Drought Index"),
          notes = c("All models include department and year fixed effects",
                    "Standard errors clustered by province are in parentheses",
                    "Models 1 -- 3: 1992-2013",
                    "Models 4 -- 6: 1992-2001"),
          out = here("Tables", "bigTable.tex"))



#####################################
##### Heterogeneity across time #####
#####################################

mRob <- plm(ln_lights ~ factor(mother)*factor(time) + factor(splinter)*factor(time) + time*(l_ln_pop + l_spei), 
            data = df1,
            model = "within",
            index=c("id_3", "year"))

# Cluster SEs
G <- length(unique(df1$id_2))
N <- length(df1$id_2)
cladjmRob <- (G/(G-1)*(N-1))/mRob$df.residual
cvcovmRob <- cladjmRob*vcovHC(mRob, type="HC1", cluster="group", adjust=T)
clmRob <- coeftest(mRob, vcov=cvcovmRob)

# Plot across time
robPlot <- data.frame(rbind(cbind(clmRob[24:44], clmRob[24:44,2], "Mother", c(1993:2013)),
                            cbind(clmRob[45:65], clmRob[45:65,2], "Splinter", c(1993:2013))))

colnames(robPlot) <- c("Estimate", "SE","Var", "Year")
robPlot[,1:2] %<>% lapply(function(x) as.numeric(as.character(x)))

# Figure 5
ggplot(robPlot, aes(y=Estimate, x=Year)) +
  geom_linerange(aes(ymin=Estimate - 1.96*SE, ymax=Estimate + 1.96*SE), lwd=1) +
  geom_vline(xintercept = 4, size=1.2, color="gray20", alpha=.5) +
  geom_pointrange(aes(ymin=Estimate - 1.68*SE, ymax=Estimate+1.68*SE), lwd=1.5,  shape=21, fill="white") +
  facet_wrap(~Var,ncol=1) + labs(x="") +
  theme_minimal() +
  geom_hline(yintercept=0, size=1.2, color="darkred", alpha=.5) +
  scale_fill_manual(values=c("white")) +
  theme(strip.text = element_text(size=11),
        panel.grid.minor.x  = element_blank(),
        panel.grid.major.x  = element_blank())
ggsave(file="robModsplot.png", path = here("Figures"), device = "png", width=7.5, height=4.5, units="in", dpi = 500)



####################################
##### Iteratively subset panel #####
####################################

# 1992-2013

df1$id <- as.numeric(factor(df1$name_3))
outMods <- vector("list", 351)
outCladj <- vector("list", 351)
clOut <- vector("list", 351)


for(i in 1:length(unique(df1$id))){
  
  outMods[[i]] <- plm(ln_lights ~ motherpost97 + splinterpost97 + time*(l_ln_pop + l_spei),
                      data = df1[df1$id != i, ],
                      model = "within",
                      effect="twoways",
                      index=c("id_3", "year"))
  
  G <- length(unique(df1$id_2[df1$id != i]))
  N <- length(df1$id_2[df1$id != i])
  outCladj[[i]] <- (G/(G-1)*(N-1))/outMods[[i]]$df.residual
  cvcov[[i]] <- outCladj[[i]]*vcovHC(outMods[[i]], type="HC1", cluster="group", adjust=T)
  clOut[[i]] <- coeftest(outMods[[i]], vcov = cvcov[[i]])
  
}

plotterM <- data.frame(matrix(nrow = 349, ncol = 2))
for(i in 1:349){
  
  plotterM[i,] <- clOut[[i]][1,1:2]
  
}

plotterM <- plotterM %>% arrange(X1) %>% mutate(lab = rownames(plotterM))

plotterS <- data.frame(matrix(nrow = 349, ncol = 2))
for(i in 1:349){
  
  plotterS[i,] <- clOut[[i]][2,1:2]
  
}

plotterS <- plotterS %>% arrange(X1) %>% mutate(lab = rownames(plotterS))


df2$id <- as.numeric(factor(df2$name_3))
outMods2 <- vector("list", 351)
outCladj2 <- vector("list", 351)
clOut2 <- vector("list", 351)

for(i in 1:length(unique(df2$id))){
  
  outMods2[[i]] <- plm(ln_lights ~ motherpost97 + splinterpost97 + time*(l_ln_pop + l_spei),
                       data = df2[df2$id != i, ],
                       model = "within",
                       effect="twoways",
                       index=c("id_3", "year"))
  
  G <- length(unique(df2$id_2[df2$id != i]))
  N <- length(df2$id_2[df2$id != i])
  outCladj2[[i]] <- (G/(G-1)*(N-1))/outMods2[[i]]$df.residual
  cvcov[[i]] <- outCladj2[[i]]*vcovHC(outMods2[[i]], type="HC1", cluster="group", adjust=T)
  clOut2[[i]] <- coeftest(outMods2[[i]], vcov = cvcov[[i]])
  
}


plotterM2 <- data.frame(matrix(nrow = 349, ncol = 2))
for(i in 1:349){
  
  plotterM2[i,] <- clOut2[[i]][1,1:2]
  
}


plotterS2 <- data.frame(matrix(nrow = 349, ncol = 2))
for(i in 1:349){
  
  plotterS2[i,] <- clOut2[[i]][2,1:2]
  
}

plotterS2 <- plotterS2 %>% arrange(X1) %>% mutate(lab = rownames(plotterS2))


# Figure 6
d1 <- ggplot(plotterS, aes(x = X1)) + 
  geom_density(alpha = .5, fill = "gray75") +
  theme_minimal() +
  geom_segment(aes(x = -0.3573849, y = -2, xend = -0.3573849, yend = 4), size = 1) +
  geom_density(data = plotterM, aes(x = X1), alpha = .5, fill = "gray75") +
  geom_segment(aes(x = -0.2327722, y = -2, xend = -0.2327722, yend = 4), size = 1) +
  theme(panel.grid = element_line()) +
  labs(x = "Estimate", title = "1992 - 2013") +
  annotate("text", x = -.379, y = 100, label = "Splinter",alpha = .5) +
  annotate("text", x = -.215, y = 100, label = "Mother",  alpha = .5)


d2 <- ggplot(plotterS2, aes(x = X1)) + 
  geom_density(alpha = .5, fill = "gray75") +
  theme_minimal() +
  geom_segment(aes(x = -0.217406, y = -2, xend = -0.217406, yend = 4), size = 1) +
  geom_density(data = plotterM2, aes(x = X1), alpha = .5, fill = "gray75") +
  geom_segment(aes(x = -0.069983, y = -2, xend = -0.069983, yend = 4), size = 1) +
  theme(panel.grid = element_line()) +
  labs(x = "Estimate", title = "1992 - 2001") +
  annotate("text", x = -0.245, y = 150, label = "Splinter",alpha = .5) +
  annotate("text", x = -0.05, y = 150, label = "Mother",  alpha = .5) +
  scale_x_continuous(breaks = seq(0,-3, -0.05))


d3 <- ggplot(plotterS, aes(x = X2)) + 
  geom_density(alpha = .5, fill = "gray75") +
  theme_minimal() +
  geom_segment(aes(x = .147, y = -2, xend = .147, yend = 40), size = 1) +
  geom_density(data = plotterM, aes(x = X2), alpha = .5, fill = "gray75") +
  geom_segment(aes(x = .143, y = -2, xend = .143, yend = 40), size = 1) +
  theme(panel.grid = element_line()) +
  labs(x = "Standard Error", title = "1992 - 2013") +
  annotate("text", x = .147 -.0006, y = 750, label = "Splinter",alpha = .5) +
  annotate("text", x = .143 -.0006, y = 750, label = "Mother",  alpha = .5) 

d4 <- ggplot(plotterS2, aes(x = X2)) + 
  geom_density(alpha = .5, fill = "gray75") +
  theme_minimal() +
  geom_segment(aes(x = .114, y = -2, xend = .114, yend = 40), size = 1) +
  geom_density(data = plotterM2, aes(x = X2), alpha = .5, fill = "gray75") +
  geom_segment(aes(x = .122, y = -2, xend = .122, yend = 40), size = 1) +
  theme(panel.grid = element_line()) +
  labs(x = "Standard Error", title = "1992 - 2001") +
  annotate("text", x = .114 - .0025, y = 750, label = "Splinter",alpha = .5) +
  annotate("text", x = .122 - .0025, y = 750, label = "Mother",  alpha = .5) 

ggsave(file="itdens.png", path = here("Figures"), device = "png", width=10, height=10, units="in", dpi = 500, arrangeGrob(d1,d2,d3,d4))



# Transfers ---------------------------------------------------------------

tran <- read.csv(here("Data", "transfersFinal.csv"))

# Figure 2

tsp1 <- ggplot(data = tran, aes(x = year, y = total / 1000000, group = province)) +
  geom_line(size = .5, alpha = .5, linetype = 2) +
  geom_line(data = tran, aes(x = year, y = mtotal / 1000000), color = "dodgerblue", size = 1.25) +
  theme_bw(base_size = 10) +
  facet_wrap(~fragcode1) +
  theme(strip.text = element_text(size=11),
        panel.grid = element_blank()) +
  labs(x = "", y = "CFA (millions)", title = "Total transfers from central government") +
  scale_y_continuous(labels = scales::comma_format()) 

tsp2 <- ggplot(data = filter(tran, province != "KADIOGO"), aes(x = year, y = (total / 1000000) / pop05, group = province)) +
  geom_line(size = .5, alpha = .5, linetype = 2) +
  geom_line(data = tran, aes(x = year, y = mtotalcap / 1000000), color = "dodgerblue", size = 1.25) +
  theme_bw(base_size = 10) +
  facet_wrap(~fragcode1) +
  theme(strip.text = element_text(size=11),
        panel.grid = element_blank()) +
  labs(x = "", y = "CFA (millions)", title = "Per capita transfers from central government",
       caption = "Dashed lines are for each individual department \n Solid lines denote mean within each fragmentation-type-year") +
  scale_y_continuous(labels = scales::comma_format()) 

ggsave(file = "transfers_time.png", 
       path=here("Figures"), device = "png", width=7, height=5.25, units = "in", dpi = 800,
       arrangeGrob(tsp1, tsp2, ncol = 1))



# DHS regressions ---------------------------------------------------------

# Infant and child mortality
dhsLights1 <- read.csv(here("Data", "dhsLightsFinal1.csv"))
# Electricity and literacy
dhsLights2 <- read.csv(here("Data", "dhsLightsFinal2.csv"))


inf_lm <- lm(inf_death ~ ln_lights + factor(year), data = dhsLights1)
child_lm <- lm(child_death ~ ln_lights + factor(year), data = dhsLights1)
elec_lm <- lm(m_elec ~ ln_lights + factor(year), data = dhsLights2)
lit_lm <- lm(reads ~ ln_lights + factor(year), data = dhsLights2)

pg_plotter <- rbind(cbind(coeftest(inf_lm)[2,][1], coeftest(inf_lm)[2,][2], "Infant Mortality"),
                    cbind(coeftest(child_lm)[2,][1], coeftest(child_lm)[2,][2], "Child Mortality"),
                    cbind(coeftest(elec_lm)[2,][1], coeftest(elec_lm)[2,][2], "Household Electrification"),
                    cbind(coeftest(lit_lm)[2,][1], coeftest(lit_lm)[2,][2], "Literacy")) %>% 
  as_tibble() %>% 
  mutate(frame = c(-0.001045773, -0.008397686,  0.085407839, 0.048721462))

# Figure 3

pgp1 <- ggplot(data = pg_plotter[1:2,], aes(y = as.numeric(V1), x = factor(V3))) +
  geom_pointrange(aes(x = factor(V3), y = as.numeric(V1), 
                      ymin = as.numeric(V1) - 1.96*as.numeric(V2),
                      ymax = as.numeric(V1) + 1.96*as.numeric(V2)), lwd=1, fill="white", shape=21) +
  theme_minimal() +
  geom_hline(yintercept=0, size=1.2, color="darkred", alpha=.5) + 
  theme(strip.text = element_text(size=11),
        panel.grid.major.x = element_blank()) +
  labs(x = "", y = "Estimate") +
  ggthemes::geom_rangeframe(aes(y = frame, x = NULL))

pgp2 <- ggplot(data = pg_plotter[3:4,], aes(y = as.numeric(V1), x = factor(V3))) +
  geom_pointrange(aes(x = factor(V3), y = as.numeric(V1), 
                      ymin = as.numeric(V1) - 1.96*as.numeric(V2),
                      ymax = as.numeric(V1) + 1.96*as.numeric(V2)), lwd=1, fill="white", shape=21) +
  theme_minimal() +
  theme(strip.text = element_text(size=11),
        panel.grid.major.x =  element_blank()) +
  labs(x = "", y = "Estimate") +
  ggthemes::geom_rangeframe(aes(y = frame, x = NULL))


ggsave(file = "dhsEsts.png", path = here("Figures"), device = "png", width=6, height=6, units="in", dpi = 500, arrangeGrob(pgp2, pgp1))





# Ethnicity ridge plots ---------------------------------------------------

# Deparment level shapefile
adm3 <- shapefile((here("Data","BFA_adm_shp", "BFA_adm3")))

# Load ethnicity rasters from side data
bobo <- raster(read.asciigrid(here("Data", "side", "bobo93.asc"), as.image = T))
dioula <- raster(read.asciigrid(here("Data", "side", "dioula93.asc"), as.image = T))
lobi <- raster(read.asciigrid(here("Data", "side", "lobi93.asc"), as.image = T))
gourmantche <- raster(read.asciigrid(here("Data", "side", "gourmantche93.asc"), as.image = T))
gouroussi <- raster(read.asciigrid(here("Data", "side", "gouroussi93.asc"), as.image = T))
mossi <- raster(read.asciigrid(here("Data", "side", "mossi93.asc"), as.image = T))
peul <- raster(read.asciigrid(here("Data", "side", "peul93.asc"), as.image = T))
senoufo <- raster(read.asciigrid(here("Data", "side", "senoufo93.asc"), as.image = T))
touareg <- raster(read.asciigrid(here("Data", "side", "touareg93.asc"), as.image = T))

# Population sums within departments
e.bobo <- extract(bobo, adm3, fun = sum, na.rm = T, df = T)
e.dioula <- extract(dioula, adm3, fun = sum, na.rm = T, df = T)
e.lobi <- extract(lobi, adm3, fun = sum, na.rm = T, df = T)
e.gourmantche <- extract(gourmantche, adm3, fun = sum, na.rm = T, df = T)
e.gouroussi <- extract(gouroussi, adm3, fun = sum, na.rm = T, df = T)
e.mossi <- extract(mossi, adm3, fun = sum, na.rm = T, df = T)
e.peul <- extract(peul, adm3, fun = sum, na.rm = T, df = T)
e.senoufo <- extract(senoufo, adm3, fun = sum, na.rm = T, df = T)
e.touareg <- extract(touareg, adm3, fun = sum, na.rm = T, df = T)

colnames(e.bobo)[2] <- "bobo"
colnames(e.dioula)[2] <- "dioula"
colnames(e.lobi)[2] <- "lobi"
colnames(e.gourmantche)[2] <- "gourmantche"
colnames(e.gouroussi)[2] <- "gouroussi"
colnames(e.mossi)[2] <- "mossi"
colnames(e.peul)[2] <- "peul"
colnames(e.senoufo)[2] <- "senoufo"
colnames(e.touareg)[2] <- "touareg"

e.bind <- e.bobo %>% full_join(e.dioula, by = "ID") %>%
  full_join(e.lobi, by = "ID") %>%
  full_join(e.gourmantche, by = "ID") %>%
  full_join(e.gouroussi, by = "ID") %>%
  full_join(e.mossi, by = "ID") %>%
  full_join(e.peul, by = "ID") %>%
  full_join(e.senoufo, by = "ID") %>%
  full_join(e.touareg, by = "ID")

e.bind$rsum <- rowSums(e.bind)
e.bind <- e.bind %>% mutate_at(.funs = funs(prop = ./rsum), .vars = vars(bobo:touareg))

adm3$ID <- as.numeric(rownames(adm3@data)) + 1
adm3eth <- full_join(e.bind, adm3@data, by = "ID")
adm3eth$name_3 <- adm3$NAME_3
adm3eth <- df1 %>% filter(year == 1993) %>% full_join(adm3eth, by = "name_3")


# Function to loop ridgeplot through all ethnic groups
rp <- function(data, eth, title){ggplot(data = data, aes_string(x = eth)) +
    geom_density_ridges(alpha = .5, aes(y = factor(split_type))) + theme_minimal() +
    labs(title = title, x = "", y = "") +
    scale_y_discrete(labels = c("No split", "Mother", "Splinter"))}

# Figure A4


ggsave(file = "ethridges.png", path=here("Figures"), device = "png", width=10, height=10, units = "in", dpi = 500,
       gridExtra::arrangeGrob(rp(data = adm3eth, eth = "bobo_prop", title = "Bobo"),
                              rp(data = adm3eth, eth = "dioula_prop", title = "Dioula"),
                              rp(data = adm3eth, eth = "lobi_prop", title = "Lobi"),
                              rp(data = adm3eth, eth = "gourmantche_prop", title = "Gourmantche"),
                              rp(data = adm3eth, eth = "gouroussi_prop", title = "Gouroussi"),
                              rp(data = adm3eth, eth = "mossi_prop", title = "Mossi"),
                              rp(data = adm3eth, eth = "peul_prop", title = "Peul"),
                              rp(data = adm3eth, eth = "senoufo_prop", title = "Senoufo"),
                              rp(data = adm3eth, eth = "touareg_prop", title = "Touareg")))

